clear all; clc; close

load EcoliHyperShockEpi.mat
i = 6;
     figure
     hold on
     
     xlabel('Time (s)')
     ylabel('L/L0')
     title(experiment{i})
     
     
     [numCells,~]=size(cell_length{i}); % numCells is number of cells
%      R = 0.000555556; % The average unperturbed growth rate
%      deltaT = R.*(t_vec{i});
%      L0 = exp(deltaT);
   
     [~,numPoints]= size(t_vec{i});
     
     locationOfShock = 32;
     
     %numPoints = locationOfShock; % when plotting only before the jump. 
 
     
    %% calculating the average length as function of time
     avg_length = zeros(numPoints,1);
     N = zeros(numPoints,1);
     for j=1:numCells
         if (~isnan(cell_length{i}(j,locationOfShock))) && (~(cell_length{i}(j,locationOfShock)==0))
           
           normalized_length = cell_length{i}(j,:)/cell_length{i}(j,locationOfShock);
           
           plot(t_vec{i},cell_length{i}(j, :)/cell_length{i}(j, locationOfShock))
             
%          normalized_length = cell_length{i}(j,:)/cell_length{i}(j,1);
%          normalized_length = normalized_length./L0 - 1;
           normalized_length(isnan(normalized_length)) = 10;
           

        %plot(deltaT,normalized_length)
        for k=1:numPoints
            if normalized_length(k) ~= 10
           
         avg_length(k) = (normalized_length(k) +  avg_length(k));
         N(k) = N(k) + 1;
            end
        end
         end
       % plot(deltaT,normalized_length)
     end
     %avg_length
     avg_length = avg_length ./ N;
    % n
     %N
     
     figure
     hold on
     
     plot(t_vec{i}, avg_length)
     xlabel('Time (s)')
     ylabel('L/L0')
     title(experiment{i})
     
     fig2pretty
    

     
        %% calculating the standard deviation corresponding to the mean above
     std_length = zeros(numPoints,1);
     N = zeros(numPoints,1);
     for j=1:numCells
         if (~isnan(cell_length{i}(j,locationOfShock))) && (~(cell_length{i}(j,locationOfShock)==0))
           
           normalized_length = cell_length{i}(j,:)/cell_length{i}(j,locationOfShock);
           
%          normalized_length = cell_length{i}(j,:)/cell_length{i}(j,1);
%          normalized_length = normalized_length./L0 - 1;
           normalized_length(isnan(normalized_length)) = 10;
           

        %plot(deltaT,normalized_length)
        for k=1:numPoints
            if normalized_length(k) ~= 10
           
         std_length(k) = std_length(k) +  (normalized_length(k) -  avg_length(k))^2;
         N(k) = N(k) + 1;
            end
        end
         end
       % plot(deltaT,normalized_length)
     end
     %avg_length
     std_length = sqrt(std_length ./ N);
    % n
     %N
     
     figure
     hold on
     
     plot(t_vec{i}, std_length)
     xlabel('Time (s)')
     ylabel('std')
     title(experiment{i})
     
     fig2pretty
     
     
     
     figure
     hold on
     
     plot(t_vec{i}, avg_length, 'b-')
     
     plot(t_vec{i}, avg_length - 2*std_length, 'r-')
     
     plot(t_vec{i}, avg_length + 2*std_length, 'r-')
     
     xlabel('Time (s)')
     ylabel('L/L0')
     title(experiment{i})
     
     fig2pretty
    
     
     %% Export arrays
       
     experiment2 = experiment;
     
     experiment2{1}(3:4) = 'to';
     experiment2{2}(4:5) =  'to';
     experiment2{3}(3:4) =  'to';
     experiment2{4}(3:4) = 'to';
     experiment2{5}(3:4) =  'to';
     experiment2{6}(3:4) =  'to';
     
     writematrix(t_vec{i},['time_'  experiment2{i} '.csv']) 
     writematrix(avg_length,['avgLength_' experiment2{i} '.csv']) 
     writematrix(std_length,['stdLength_' experiment2{i} '.csv']) 
     
 %%
%      fullSamples = [];
%      for j = 1:numCells
%          
%        if  (~isnan(cell_length{i}(j,locationOfShock))) && (~(cell_length{i}(j,locationOfShock)==0))
%             fullSamples = [fullSamples j];
%        end
%      end      
%      
%      figure 
%      hold on
%      
%      
%      for j=fullSamples
%         
%      end
%    